Liquid Polymorphism and Density Anomaly in a Lattice Gas Model 
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We present a simple model for an associating liquid in which polymorphism and density anomaly 
are connected. Our model combines a two dimensional lattice gas with particles interacting through a 
soft core potential and orientational degrees of freedom represented through thermal "ice variables". 
The competition between the directional attractive forces and the soft core potential leads to a phase 
diagram in which two liquid phases and a density anomaly are present. The coexistence line between 
the low density liquid and the high density liquid has a positive slope contradicting the surmise that 
the presence of a density anomaly implies that the high density liquid is more entropic than the low 
density liquid. 
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Water is one of the most mysterious materials in na- 
ture. It exhibits a number of thermodynamic and dy- 
namic anomalous properties [l|-|3|], such as the maxi- 
mum as a function of temperature both in density and 
in isothermal compressibility in the liquid phase. It 
has been proposed some time ago that these anomalies 
might be associated with a critical point at the terminus 
of a liquid-liquid line, in the unstable supercooled liq- 
uid region^, at high pressures, following the suggestion, 
based on varied experimental data |5|, of a thermody- 
namic singularity in supercooled water, around 228K and 
at atmospheric pressure. Inspite of the limit of 235K be- 
low which water cannot be found in the liquid phase with- 
out crystallization, two amorphous phases were observed 
at much lower temperatures [||. There is evidence, al- 
though yet under test, that these two amorphous phases 
are related to fluid water @ • 

Notwithstanding of its confirmation for metastable wa- 
ter, interest in liquid polymorphism arose, and the coex- 
istence of two liquid phases was uncovered as a possibility 
for a few other both associating and non-associating liq- 
uids. Notable exam ples include liquid metals silica 
[Tof. phosphorus [ITI||l2j and graphite 0. The relation 
between liquid polymorphism and density anomaly has 
been a subject of debate in recent theoretical literature 
0. 

From a microscopic point of view, water anomalies 
have been interpreted qualitatively, since|l5|. in terms 
of the the presence of an extensive hy drog en bond net- 
work which persists in the fluid phase [lg. The hydro- 
gen bond net deformation under temperature was repre- 
sented in a number of minimal models for water. The 
main strategy has been to associate the hydrogen bond 



disorder with bond [HUH or site [SJ Potts states. 
In the former case coexistence between two liquid phases 
may follow from the presence of an order-disorder tran- 
sition and a density anomaly is introduced ad hoc by 
the addition to the free energy of a volume term propor- 
tional to a Potts order parameter. In the second case, 
it may arise from the competition between occupational 
and Potts variables introduced through a depedency of 
bond strength on local density states. 

We propose a description also based on occupational 
and orientational degrees of freedom. However, for the 
orientational part we employ a modification of the ther- 
mal version pHp^| of the ice modeljl^], so successful in 
the description of ice entropy. Competition between the 
filling up of the lattice and the formation of an open four- 
bonded orientational structure is naturally introduced in 
terms of the ice bonding variables and no ad hoc intro- 
duction of density or bond strength variations is needed. 
We thus consider a lattice gas on a triangular lattice with 
sites which may be full or empty. Besides the occupa- 
tional variables, Cj, associated to each particle i there 
are six other variables, -n , pointing to neighboring sites 
j: four are the usual ice bonding arms, two donor, with 
r? J = 1, and two acceptor, with = — 1, while two ad- 
ditional opposite arms are taken as inert (non-bonding), 
t* 3 = 0, as illustrated in Fig. 1. Therefore each occupied 
site is allowed to be in one of eighteen possible states. 

Two kinds of interactions are considered: isotropic 
"van der Waals" and orientational hydrogen bonding. An 
energy — v is attributed to each pair of occupied neighbor- 
ing sites that form a hydrogen bond, while non-bonding 
pairs have an energy, — v + 2u (for u > 0). The overall 
model energy is given by 



E = ^ o-iCTj{-v + u[2 
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where Ui = 0, 1 are occupation variables and r l y 
represent the arm states described above. 
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FIG. 1: The model orientational state: four bonding (donor 
and receptor) and two non-bonding arms 







FIG. 2: Low density liquid, LDL with density 3/4 (top), 
and high density liquid, HDL witn density 1 (botton) on the 
triangular lattice. 



Inspection of the model properties allows the predic- 
tion of two ordered states, as shown in Fig. 2. One 
of them has lower density ,p = 0.75, and energy den- 
sity given by e = E/N = — 3v/2. The other state has 
higher density ,p — 1, and energy density e = — 3v + 2u. 
The addition of an external chemical potential p may 
favor one or the other of the two ordered states. At 
zero temperature, the low density liquid (LDL) coexists 
with the high density liquid (HDL) at chemical potential 
fj,/v = —6 + 8u/v, obtained by equating the grand po- 
tential density (or pressure) associated with each one of 
these phases. Similarly the coexistence pressure at the 
zero temperature is given by p/v = —3 + 6u/v. Besides 
these two liquid states, a gas phase is also found and it 
coexists with the low density liquid at chemical potential 
fj,/v = —2 and pressure p = 0. The condition for the 
presence of the two liquid phases is therefore u/v > 0.5. 

The model properties for finite temperatures were ob- 
tained through Monte Carlo simulations in the grand- 
canonical ensemble using the Metropolis algorithm. Par- 
ticle insertion and exclusion were tested with transition 
probabilities given by 



w(insertion) 
w(exclusion) 

w(insertion) 
w(exclusion) 



exp(—A(f>) and 
1 if A0>O 
or 

1 and 

exp(+Acj)) if AcfxO 




(2) 



FIG. 3: Pressure vs. density isotherms for differ- 
ent temperatures. 



with A(j> = exp{/3(e partic ; e -^)-ln(18)} where e par ticie is 
the energy of the particle included. Since the empty and 
full sites are visited randomly, the factor 18 is required 
in order to guarantee detailed balance. 

Simulational data were generated both from fixed tem- 
perature and fixed chemical potential simulations. Some 
test runs were done for L=4, 10 and 20. A detailed study 
of the model properties and the full phase diagrams was 
undertaken for an L=10 lattice, at fixed u/v = l[24|. 
Runs were of the order of 10 6 Monte Carlo steps. 

The three phases obtained at zero temperature are 
present for low temperatures, as can be seen in the 
isotherms of Figure 3. The model exhibits two first order 
phase transition lines, gas-LDL and LDL-HDL, respec- 
tively. 

In order to obtain the complete phase diagram, in- 
cluding the two critical points, and to check for density 
anomalies, pressure was computed by numerical integra- 
tion of the Gibbs Duhem equation, SdT — VdP + Nd[i — 
0, at fixed temperature. Integration was carried out from 
effective zero density, at which pressure is zero, to obtain 
P(p, T) isotherms. 

The pressure isotherms show that an inversion of the 
behavior of density as a function of temperature takes 
place at intermediate pressures, in the LDL phase. At 
smaller pressures ,p ~ 1, density decreases with temper- 
ature, whereas at higher pressures ,p ~ 3, density in- 
creases with temperature. This yields a density anomaly 
in the higher range of pressures, which we illustrate in 
Fig 4. 

Finally, from a large set of temperatures, we build 
up the pressure versus temperature coexistence curves 
shown in Fig. 5. The line of maximum densities is also 
shown. The liquid-liquid coexistence line has a positive 
inclination, except at very low temperatures (the zero 
temperature points do not come from simulations, but 
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FIG. 4: Density Anomaly for different pressures 
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FIG. 5: Phase Diagram. The open circles repre- 
sent the temperature of maximum density, TMD, 
line, the filled circles represent the LDL-HDL co- 
existence that ends at the critical point C and 
the filled squares represent the gas-LDL coexis- 
tence line ending at the critical point C. 



from equating the enthalpy densities). From the Clapey- 
ron condition, ^} CO ex — j the positive slope implies 
that within our model the HDL phase has lower entropy 
than the LDL phase. 

Inside each phase, the density anomaly can be related 
to the behavior of entropy as a function of pressure. From 
thermodynamics, a negative thermal expansion coeffi- 
cient a = (^)pimplies a positive gradient of entropy 
with pressure, since (§f) P — ~(§§)t- This property has 
been thought [l4| to imply that the presence of a density 
anomaly would lead to a high entropy high density phase, 
and therefore to a negative slope of the coexistence line, 



as is true for the ice fusion line. The present model proves 
that this assumption is misfounded and that this is not 
a general behavior. 

What we have here is the following: on the low den- 
sity side, the thermal expansion coefficient is negative, 
whereas on the high density phase it is positive, as can 
be gathered from the pressure-density isotherms. The 
positive slope of the coexistence line implies, by Clausius- 
Clapeyron, that the high density phase is the lower en- 
tropy phase. Thus, at constant temperature, entropy 
increases with pressure up to the coexistence line, drops 
discontinuously across this line, and then decreases with 
pressure, as in any normal liquid. Therefore the sign 
of entropy variation across the coexistence line may be 
either positive, as in this model, or negative, as in the 
fusion of ice, following, in both cases, the high pressure a 
sign. 

The model proposed is a truly statistical model which 
includes orientational and occupational variables, and 
guarantees the local distribution of hydrogens on molec- 
ular bonds, without the need of increasing the volume 
artificially or introducing artificial orientational vari- 
ables. Inspite of the absence of an orientational order- 
disorder transition |2^ . the model presents liquid- liquid 
coexistence, with positive inclination in the pressure- 
temperature plane, accompanied by a line of maximum 
density, on the low density side, a feature expected for 
water. Besides, this study points out to the fact that the 
presence of a density anomaly, with a < 0, on the low 
temperature side, and as a consequence, (§|)t > 0, does 
not imply a negative slope of the liquid-liquid line, con- 
trasting with the results for most studies of metastable 
liquid-liquid coexistence in models for water, which sug- 
gest a transition line with negative gradient [25| . 

The presence of both a density anomaly and two liquid 
phases in our model begs the question of which features 
of this potential are responsible for such behaviour. Av- 
eraged over orientational degrees of freedom, our model 
can be seen as some kind of shoulder potential, with 
the liquid-liquid coexistence line being present only for 
a repulsive, soft-core potential. The same was indeed 
observed for continuous step pair potentials 14] 26], for 
which, however, the density anomaly is absent. On the 
other hand, a density anomaly seems to be associated 
with smooth soft core potentials [27] 28], which would 
be hidden, in our model, in the orientational degrees of 
freedom. 

In summary, we have found that a lattice gas with 
orientational ice-like degrees of freedom can generate a 
density anomaly and a liquid-liquid phase boundary with 
positive slope. 
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